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Abstract 

As well known, Boltzmann-Gibbs statistics is the correct way of thermostatistically approaching 
ergodic systems. On the other hand, nontrivial ergodicity breakdown and strong correlations 
typically drag the system into out-of-equilibrium states where Boltzmann-Gibbs statistics fails. For 
a wide class of such systems, it has been shown in recent years that the correct approach is to use 
Tsallis statistics instead. Here we show how the dynamics of the paradigmatic conservative (area¬ 
preserving) standard map exhibits, in an exceptionally clear manner, the crossing from one statistics 
to the other. Our results unambiguously illustrate the domains of validity of both Boltzmann-Gibbs 
and Tsallis statistics. 

PACS numbers; 05.20.-y , 05.10.-a , 05.45.-a 


*Electronic address; ugur.tirnakli@ege.edu.tr 
^Electronic address; ernesto@ufba.br 


1 



I. INTRODUCTION 


Exponential and Gaussian distributions are signatures of the Boltzmann-Gibbs statistical 
mechanics. These distributions are those that maximise the Boltzmann-Gibbs entropy and 
ensure the equilibrium state. The Maxwell distribution is an instance of the equilibrium dis¬ 
tribution for the velocities of molecules in an ideal gas. The underlying mathematical reason 
for this is the existence of the standard Central Limit Theorem (CLT) [I] . On the other hand, 
due to ergodicity breaking, some systems remain indefinitely trapped into non-exponential 
and non-Gaussian distributions, and thus achieve out-of-equilibrium quasi-stationary states. 
The g-exponential and the g-Gaussian distributions are functions associated with some of 
these quasi-stationary states and they are the maximising distributions for the non-additive 
Tsallis entropy given by = A:(l — ’Yhi'Pi) / (? ~ 1 ) [ 2 ]- This feature permits to describe 
these special non-equilibrium states with the same formal framework of the equilibrium 
thermostatistics, known as Tsallis statistics [3], and this general picture is reduced to the 
equilibrium one if the parameter g attains a special limiting value (g ^ 1). In this case, 
the underlying mathematical mechanism is the generalized CLT illH], which states that the 
stable limit distributions for a certain class of systems in such quasi-stationary states are 
g-Gaussians. Therefore, the role of g-Gaussians in Tsallis statistics is basically the same as 
that of Gaussians in Boltzmann-Gibbs statistics. In this work we show, for the first time, 
that these two cases coexist in the classical standard map, and discuss the necessary condi¬ 
tions under which one case prevails over the other one. This neatly illustrates the respective 
domains of validity of Boltzmann-Gibbs and of Tsallis statistics. 

Non-Gaussian distributions, particularly g-Gaussians, have been observed in nature in 
several experimental, observational and model systems j3]. Impressive experimental exam¬ 
ples include (i) a high dimensional dissipative system where the probability densities of 
velocity differences measured in a Couette-Taylor experiment for a fully developed turbu¬ 
lence regime pm, (ii) transport properties of cold atoms in dissipative optical lattices Pi 
and (iii) transverse momentum spectra of hadrons at LHG experiments [10]. As observa¬ 
tional works, for small bodies in the Solar System, particularly asteroid rotation periods 
and diameters m and distribution of meteor showers [T2| can be given. At a larger scale, 
the rotation curve for the M33 Triangulum Galaxy has been successfully analyzed in the 
same sense [I3|. Among model systems, one of the paradigmatic dissipative low dimen- 
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sional model, the logistic map, has been numerically investigated and g-Gaussians have 
been found as the chaos threshold is approached using the band splitting structure obeying 
the Huberman-Rudnick scaling law [TH - fTt)] . 

g-Gaussians have also been recently observed in a conservative high dimensional model 
HZ]. In the a-XY model, i.e., a system of N classical localized planar rotators with two- 
body interactions and periodic boundary conditions, the potential is assumed to decay with 
distance as l/r“, and « > 0 is the parameter that controls the range of the interactions, 
short-range for a/d > 1, and long-range for 0 < a/d < 1 (d is the spatial dimensionality of 
the system). 

Recently a generalization of the conservative one-dimensional Fermi-Pasta-Ulam model, 
properly modified to account for linear and nonlinear long-range interactions, has been 
analyzed. The range of the interactions is controlled in the same way as for the a-XY model 
just mentioned. Ordinary Gaussians are observed when short-range interactions (a > 1) are 
present, and g-Gaussians are observed when long-range interactions (0 < a < 1) are present 
[IH]. It has been found that the maximal Lyapunov exponent A asymptotically decreases 
as in a rather similar way of that observed in [19] for the a-XY model and the 

g-Gaussian distributions that emerge are characterized by the parameter g that depends on 
a. 

All these systems appear to share in common the following scenario: ergodicity in a region 
is characterized by the largest Lyapunov exponent A and two regimes shall be distinguished 
in the thermodynamic limit (number of particles N oo). (i) Strongly chaotic regime cor¬ 
responds to a large positive Lyapunov exponent, where the system is ergodic. The dynamics 
of the system evolves to an equilibrium state described by Boltzmann-Gibbs statistical me¬ 
chanics, with exponential or Gaussian distributions (according to the considered dynamical 
variable); (ii) Weakly chaotic regime corresponds to a very small positive Lyapunov expo¬ 
nent (A ~ 0), where the system behaves for a very long time as non-ergodic. Distributions 
of the dynamical variables are not exponential or Gaussians, and the Boltzmann-Gibbs 
framework is not suitable for this case. g-Gaussian distributions have been observed for 
this case by proper time and ensemble averages [20|. These distributions are obtained by 
maximisation of the nonadditive entropy Sg PEHEg, which is a strong indication that 
these systems are connected to nonextensive statistical mechanics [Hj. They may be written 
as P{u) oc expg(—Ru^), where R > 0 is the Lagrange parameter and the g-exponential is 
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given by exp^n = [1 + (1 — with [^4]+ = max{0,74} and its inverse, g-logarithm, 

is defined by In^n = {v}~^ — 1)/(1 — g). The ordinary Gaussian, exponential and logarithm 
functions are respectively recovered in the limit g ^ 1. 

In this paper we consider the standard map, that is a paradigmatic low dimensional 
conservative (area-preserving) model, and we follow the averaging procedure originally used 
for the logistic map, as described in [I1HI6]. As will be discussed in detail below, this 
paradigmatic model offers an excellent medium for us to analyse both regimes explained 
above and to establish a connection between these regimes where the system is ergodic and 
non-ergodic. 

II. THE MODEL AND RESULTS 

The standard map is defined as [2HI El! 

Pi+i = Pi - K sin Xi ; Xj+i =Xi+ Pi+i (1) 

where p and x are taken as modulo 27r. This map has very rich properties depending on the 
map parameter K. Here, we will focus on four representative cases whose phase portraits are 
given in Fig. The two extreme cases are K — 0.2 and K = 10, one of which represents the 
domination of the phase space with the stability islands and the other is clearly an example 
of the invasion of the full phase space by the chaotic sea. On the other hand, the other two 
cases in between, namely, K = 0.6 and K = 2, are good examples in order to see how these 
regions with stability islands and chaotic sea merge in the available phase space. It is clear 
that if the system starts from an initial condition located on one of the archipelagos (given 
by the same color), it will stay forever in the same archipelago, whereas if it starts from 
somewhere in the chaotic sea, the iterates will cover the whole chaotic region. 

At this point, we need to calculate the largest Lyapunov exponent of these cases using 
the Benettin algorithm [2S] but this calculation is to be done very carefully. Generally, 
calculating the Lyapunov exponent by taking an ensemble average would not be exactly 
correct here since the contributions coming from the initial conditions of stability islands 
are much smaller than the ones coming from the chaotic sea. Therefore, making an ensemble 
average would not reflect the correct behaviour of the system. In order to reflect the correct 
behaviour, we prefer to plot the largest Lyapunov exponents as given in Fig. where we 
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FIG. 1: (Color online) Phase portrait of the standard map for 4 representative K values. In each 
case, black dots represents the region of chaotic sea in the available phase space and all other colors 
represent different stability islands. 

calculate the exponent of each initial condition separately over the whole phase space and 
the magnitude of the exponents are given by a color map. As seen in the figure, the case 
A' = 0.2 represents a Lyapunov spectrum in which all results are extremely close to zero 
(black dots), whereas the case K = 10 conversely exhibits a spectrum where all results 
are largely positive (yellowish dots). This means that, in the former case (latter case), the 
whole phase space is dominated by the stability islands (chaotic sea). On the other hand, 
the other two cases, K = 0.6 and K = 2, are good examples where the phase space consists 
of both stability islands and chaotic sea. This way of representing the Lyapunov spectrum 
allows us to see clearly the portions of the whole phase space where the system is ergodic 
and non-ergodic for a given K value. 

Now we can analyze the limit distributions of the standard map for these representative 
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FIG. 2: (Color online) Lyapunov exponent results of the phase portrait of the standard map. 
The same representative K values are used. For each case, Lyapunov exponents are calculated for 
200000 initial conditions. In the calculation, each initial condition is iterated 10^ times. 


K values. We define the variable 


y 


■■= - (x)), 


( 2 ) 


i=l 


where the average (• • •) is calculated as time average taken over not only a large nnmber of 
T iterations, but also a large number of M randomly chosen initial values, namely, 

M T 

(3) 

i=i «=i 

and calculate the probability distribution of y, namely P{y), for any given K parameter. 

Let us start with the result of the case K = 10, where the probability distribution is 
expected to be a Gaussian since for this case the phase space is totally a chaotic sea, which 
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makes the whole system ergodic. The result is given in Fig. where a clear Gaussian is 
easily seen as expected. It should also be noted that the stable limit distribution is obtained 
quickly. 
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FIG. 3: (Color online) Normalized probability distribntion function for the case K = 10 with 
T = 2^®. 

Now we can investigate the case K = 0.2, where the probability distribution is expected 
to be a non-Gaussian due to the change of the phase space from being totally chaotic to 
totally consists of stability islands, which makes the whole system non-ergodic. The result 
is given in Fig. where, instead of a Gaussian, now a clear g-Gaussian is observed with 
q = 1.935. In this case, the stable limit distribution happens to be achieved slowly but at 
the level of T = 2^^ it has already been reached. We have checked it with T = 2^® and 
verified that the distribution does not change in the displayed region. We also plot the 
same data as g-logarithm of the probability distribution in Fig. in order to see whether 
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it is a straight line or not. For three different regions (i.e., the region including the tails, 
intermediate region and the central part), the straight line is well approached. The fact that 
we observe straight lines in all scales excludes other distributions that are asymptotic power 
laws, like Levy distributions. 



FIG. 4: (Color online) Normalized probability distribntion function for the case K = 0.2 with 
T = 2^^. In the Inset, the central part is zoomed for a better visnalization. 

In order to better illustrate this tendency, we also perform another test by taking a K 
value where the stability islands and chaotic sea coexist, i.e., the case K = 2. For this case, it 
is evident from Fig. that the region of chaotic sea with large positive Lyapunov exponents 
and the region of stability islands with Lyapunov exponents close to zero can easily be 
detected. This means that the system is ergodic within some portion of the phase space, 
whereas it is indeed non-ergodic within some other portion. Therefore we can check our 
previous findings using these portions separately. If we use initial conditions all taken from 













c 


II 

a" 



[y P(0)]" 


FIG. 5: (Color online) g-logarithmic representation of the normalized probability distribntion 
(a) for the tails, (b) for the intermediate region and (c) for the central part of the case K — 0.2. 

the portion where the system is ergodic (non-ergodic), we expect to see the same distribution 
function we have found before, namely the Gaussian (g-Gaussian with q = 1.935). In fact, 
this is exactly what we see in Fig. which nicely corroborates our results given in Fig. 
and Fig. 

Finally we will be interested in another interesting question: what happens to the prob¬ 
ability distribution if we do not take the portions of the phase space separately where the 
system is ergodic and non-ergodic but consider initial conditions coming from the whole 
phase space in the calculation of the probability distribution. This is really worth analysing 
since in this case one would expect a competition between initial conditions coming from 
the regions of the available phase space where the system is ergodic and non-ergodic and 
therefore between Gaussian and g-Gaussian behaviour. Needless to say, as the region in the 
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phase space where the system is non-ergodic diminishes (like the case K = 10), Gaussian 
distribution will win, whereas the winner will be g-Gaussian as the region where the system 
is ergodic shrinks (like the case K = 0.2). We notice that, if these regions coexist together, 
then this competition between Gaussian and g-Gaussian can be modelled as 

= a expg{-/3q[yP{0)]‘^) + (1 - a) exp{-P[yP{0)f). (4) 

We check this hypothesis using our two appropriate cases, namely, K = 0.6 and K = 2. The 
results are given in Fig. where a corroboration can be seen at different scales. 


III. CONCLUSIONS 

The phase space of the standard map presents regions of positive Lyapunov exponents 
coexisting with regions of zero Lyapunov exponents. The positive Lyapunov regions present 
mixing and thus the system is ergodic in those regions. For sufficiently low values of the con¬ 
trol parameter K, the phase space is almost entirely dominated by zero Lyapunov behaviour 
and the distributions (obtained through time averaging, along the lines of central limit theo¬ 
rems) are g-Gaussians. As the value of K increases the measure of the zero Lyapunov regions 
decreases, and we see a continuous crossing, expressed by the parameter a in Eq. 0 , be¬ 
tween g-Gaussians distributions (Tsallis statistics) and Gaussian ones (Boltzmann-Gibbs 
statistics) with Pg ^ oo and /3 —>■ tt as W —>■ oo (that is equivalent to a —>■ 0). Remarkably 
enough, the distributions originated from initial conditions taken inside the region of islands, 
instead of over the entire phase space, yield one and the same value g = 1.935, independently 
on whether we consider one or many of these regions, and independently from K. Initial 
conditions taken within the chaotic sea always yield Gaussians with {3 = tt. The variance of 
g-Gaussians with 5/3 < g < 3 diverges, though they have finite width. The A^-fold convolu¬ 
tion product of independent (or quasi-independent) g-Gaussians would asymptotically yield 
Levy distributions [21]. Fig. [^neatly shows that this is not the case for the standard map: 
indeed, the actual time-averaging involves strong correlations. 
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FIG. 6: (Color online) Normalized probability distribution function for the case K = 2. In the 
calculations, all initial conditions are taken from the region of (a) chaotic sea and (b) stability 
islands. 
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FIG. 7: (Color online) Normalized probability distribution function for the cases (a) K = 0.6 
and (b) K — 2. In the calculations, initial conditions are randomly taken from the whole available 
phase space. 
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